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SUMMARY 

An extension of the approximation factorization (AF) scheme to multi-block 
grids Is studied. An Innovative method Involing zonal decomposition using 
interfacing grids is investigated. Suitably matching the solutions at a block 
boundary for an incompressible flow problem it is demonstrates that this 
method is feasible for both contiguous and noncontiguous interfacing grids. 
Remaining issues for extending this method to transonic flow problem are 
discussed. The flow solution algorithm of a fully implicit AF code (TWINGB) 
is investigated for adaptation to block-structured grids . The analysis shows 
that the TWINGB code is sensitive to grid effects. A method to improve the 
TWINGB flow solution algorithm is presented to reduce the grid dependency of 
the code. Test case studies for subsonic and transonic flow over parabolic 
arc airfoils and wings suggest that the improved method should be adapted for 
extending TWINGB more effectively to flow problem involving arbitrary grids 
with block structure. 
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I. INTRODUCTION 


The objective of the present study is to explore means for extending 
approximate factorization (AF) schemes (reference 1) to the treatment of 
transonic flows requiring use of complex, multi-block grids (reference 2) such 
as dipicted in Figure 1. It has been demonstrated on a single-block ring grid 
that the fully implicit nature of the AF scheme (reference 1) significantly 
improves the convergence speed of the flow solution relative to successive 
line overrelaxation (SLOR) scheme. It is highly desirable to extend this 
scheme to block-structured grid. The primary issues involved center about 
devising means for extending the transonic finite difference algorithms beyond 
the simplicity inherent in a simple, single-block grid. 

The main advantage of block-structured grids is their adaptability to complex 
configurations. However, special issues are introduced in the flow solution - 
algorithms such as fictitious corners and nonanalytic block boundaries. These 
algorithm issues have been studied in two-dimensions (reference 3) and 
three-dimensions (reference 2) when the flow is solved by a finite volume 
method using a SLOR scheme. 

Two innovative approaches for solving the transonic flow in a multi-block grid 
were explored. The first approach examines a method involving “zonal 
decomposition" wherein block boundaries are treated as true boundary surfaces 
separating interfacing grids. The issues investigated involve techniques for 
matching solutions at a block boundary. A feasibility study was completed and 
the results are presented. 

The second approach (Figure 2) involves overlapping grids for differencing 
across a block boundary near an artificially induced coordinate singularity 
occurring at a fictitious corner. This approach selects a set of neighboring 
nodes for the fictitious corner (node 3, see Figure 2c) such that the 
resulting physical cells (cell 1342 and cell 3564) for node 3 are 
topologically the same as any other node on the airfoil surface. 




' — Fictitious comer 

Physical space 



Computational space 


Figure 1. Block Structuring 
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(c) Near a fictitious 
corner 



Figure 2 . A Typical Block-Structured Grid. 
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To explore the second approach, an AF-based transonic code,, called TWINGB, was 
modified for the treatment of flow using multi-block grids.- During the course 
of this study. It was discovered that the flow solution algorithm In TWINGB Is 
sensitive to grid effects. A method to Improve the flow solution algorithm In 
TWINGB to reduce Its grid dependency has been derived. The modified TWINGB 
has been evaluated on test cases with a typical parabolic arc airfoil and wing 
using block-structured grids. Results of the study Indicate a further need to 
Improve the flow solution algorithm. 

The block-structured grid Is discussed in Chapter II. The fully implicit AF 
scheme Is given In Chapter III. The method of solution employing zonal 
decomposition Is presented in Chapter IV. The study of grid effects for the 
AF-based TWINGB code Is summarized In Chapter V. Results for adaptation of 
TliJINGB to block-structured grids are presented In Chapter VI. Conclusions and 
recommendations are given In Chapters VII and VIII. 
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II. BLOCK-STRUCTURED GRID 


The construction of a suitable grid system for complex 3D configurations, such 
as a wing/body /nacelle, is a necessary step for computing the corresponding 
transonic full potential flow. Two approaches have been available based on 
Thompson's (reference 4) surface-adapted coordinate concept. A limited 
approach maps the flow domain surrounding a 3D configuration into a single 
rectangular box. This approach has been successful for simple geometries 
(reference 5) but cannot be effectively applied to complex configurations with 
multiple components. The more general approach divides the computational 
domain into multiple rectangular blocks where the configuration itself is also 
represented by a set of blocks whose structure follows the nature lines of the 
configuration (Figure 1). This results in a multi-block grid system which is 
adaptable to complex configurations and can produce good grid quality near 
physical corners. A comparison of typical single block and multi-block grids 
is given in Figure 3. 

Generally, a single-block grid structure will not be suitable for 
geometrically complex boundary configurations. Thus, the flow field must be 
divided into multiple blocks to accommodate slope discontinuities of the 
boundary surfaces and provide good boundary fitting behavior (Figure 4). 

A typical block-structured grid generation (Figure 5) process is desribed 
below: 

1. Define the overall block structure according to the natural lines of 
the configuration. 

2. Generate ID grid along the block perimeters (perimeter 
discretization). 

3. Generate 2D grid covering the block surfaces using the information 
obtained in Step 2 as boundary data for the surface grid generation. 

4. Generate 3D grid covering the volume grids filling each block using 
the information obtained in Step 3 as boundary data for the volume 
grid generation. 
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Figure 3. Comparison of Grids Near an Airfoil. 
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Figure 4. Comparison of Grid Structure. 
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Thompson's method of generating a boundary-fitted curvilinear coordinate 
system requires that the curvilinear coordinates be solutions of a system of 
elliptic partial differential equations in the physical space subject to 
Dirichlet boundary conditions on all boundaries. One (curvilinear) coordinate 
is -specified to be constant, on each of the boundaries and a monotonic 
variation of the other coordinate along each boundary is specified. By 
interchanging the role of the Cartesian coordinates and the transformed 
coordinates, one arrives at a quasi-linear elliptic system for the Cartesian 
coordinates in the transformed space. 

The present approach simplifies the formulation by requiring the Cartesian 
coordinates x * (x, y, z) to be the solutions of three independent linear 
equations defined in the computational coordinates ( C , n, ^ 

where the coefficients A to F are constants or specified functions used for 
grid control. The coefficients A, B, and C provide means to locally rescale 
the computational coordinates. Along the block boundaries, these three 
coefficients may be evaluated knowing the grid distributions, e.g. 

A « (x^x x^)^ (2) 

The coefficients D, E, and F are related to the source terms in Thompson's 
nonlinear approach (reference 4). This can be seen by comparing equation (1) 
with the quasi-linear equations for 3D grid generations (more details can be 
seen in reference 2). These three coefficients can be extracted effectively 
from the perimeter discretization of the block boundary (reference 6) taking 
the limiting form of equation (1). For example, D may be extracted from the 
relation. 

=0 (3) 

The coefficients D, E, and F are linearly interpolated on a block boundary 
surface using the block perimeter value. The grid control coefficients A to F 
are then linearly interpolated inside a block using its boundary values. 
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III. AF2 SCHEME 


A fast, fully implicit AF algorithm to solve the conservative, transonic, full 
potential equation in 3D was presented by Holst (reference 1) in a 
single-block ring grid. General guildeines for the construction of implicit 
AF schemes follows from the two-level iteration procedure. 

NC" + uLd" « 0 (4) 

where «* d” ^ - d” is the correction, Ld” is the residual of the nth 
level velocity potential solution (d*^), and w is a relaxation parameter. 

In the AF approach, N is factorized by several factors, usually, two factors 
for 2D algorithms and three factors for 3D algorithms. These factors are 
chosen so that their product is an approximation to L, only simple matrix 
operations are required to ensure computational efficiency, and the overall 
scheme is stable. 

Stability in the supersonic regions of flow has been achieved by adding an 
artificial viscosity term using an upwind bias of the density coefficient. 


A. Governing Equations 

The 3D full-potential equation in strong conservation-law form is given by 


p - [1 - ^ Wi 4^)3 (5b) 

The density (p) is normalized by the stagnation density (p^), the three 

velocity components (d^* dy^ and d^) are normalized by the critical 

sound speed; x, y, and z are Cartesian coordinates in the streamwise, spanwise. 
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and vertical directions respectively, and y is the ratio of specific heats. 

Equation (5) is transformed from the physical space (x, y, z) to the 
computational space ( n, S ) using a general independent variable 
transformation. This transformation, denoted by 

5 “ 5 (x,y,z) n » n(x,y,z) 5 » 5 (x,y,z) ( 6 ) 

maintains the strong conservation-law form of Equation (5). The full 
potential equation written in the computational domain is given by 

(pU/D)j + (pV/D)^ + (pW/D)^ « 0 (7a) 


D ® ( VnxVS) “ det(J) 


(7b) 


q 


2 



•u' 

1 

'V 

a i 

V 


4 > 




n 


..w . 




(7c) 


p 


[1 


(7d) 


U = 

V = Vn-v^ ( 8 a) 

W = V5*Vd 
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and 



U (V or W) is the contravariant velocity component normal to the %(n or 5 ) 
plane, to Ag are matrix quantities, and J is the Jacobian matrix of the 
transformation. The expression in Equation (8) can be easily evaluated using 
the relation 



Equation (9) is valid so long as the Jacobian matrix is nonsingular. 

Equations (7) and (8) can be simplified if all constant n surfaces coincide 
with constant y planes. That is; 


y » 0 = 0 (10) 

The metric quantities of Equations (8) and (9) become: 
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, - 5^ * 5^ * ■s' 

3 ^ X ^2 


55 + 3r5+j< 

^x -^x y -'y •'z ^z 


A 4 = Sy ny 

*6 = Sy "y 


(11a) 


D . l/(Xj Zj - Xj Zj )y^ 


and 


*x - “(y, H • 


^x " -“(I'r 


Jy . D(x^ z^ - x^Zj ), 


^Z = -“<*5 1'4> 


"x ' “ "y - 


5y- -*sV 


5j= “("sV 


"z = ° 


(12) 


B. 


Flow Solution Algorithm 


Let i,j,k be the subscripts indicating position in the 5 , t), 5 computational 
mesh. Referring to a point i,j,k, a second order accurate finite difference 
approximation to equation (7a) is given by: 




J.k+I" ° 


The operator (a or «-) is first order accurate backward difference 
> n > 

operator in ^ (n or- 5 ) direction. These operators are defined by: 


(13) 



^',j»k 

» ( 

^i.j.k " ^ 

^i-l,j,k 


^’J,k 

» ( 

^i.j.k ■ ^ 

^,j-l,k 


^,j.k 

= ( 

^i.j,k ■ ^ 

^".j,k-l 


( 14 ) 
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It should be noted that equation (14) provides the divided difference formula 
since the standard divider » Aji» and a^are all equal to one and have been 
eliminated. The density is calculated using equation (7d). Values of density 
are computed. and stored at half-I point in the computational mesh. The n» 
and ^ derivatives of required for computing the density at the half-I point 
i+l/2,j,k, are approximated by: 


N+l/2.j,k 

d 

’’i+l/2.j,k 

^i+l/2.j.k 


^i+l,j,k - ^i,j,k 


°*25(^i+l,j,k -^i+lj-l,k 


^i,j+l,k - ^i J-l,k^ 


0.25(di+i j^k+l - ^i+lJ-l.k-1 


^i,j+l,k+l - ^i,j-l,k-l ^ 


(15) 


The contravariant velocity components, j ic» j+l/2 k* 

Wi j k+i/2 computed by standard second order accurate finite difference 
formulas, e.g.. 


^i+l/2,j,k " \+i/2,j,k^^i+l/2,j,k 

^i+l/2,j,k ’’i+l/2,j,k 
+ A «L 

^i+l/2,j,k^i+l/2,j,k (16) 

where the approximation formulas in equation (15) are used to provide the 
difference approximations in 5 , n and 5 for i. 
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In order to maintain stable supersonic region, the following artificial 
viscosity term was added to equation (13) 

-A^(vp - An(yp (17) 

where v is an artificial viscosity coefficient to be defined later. The 
addition of such an artificial viscosity term is equivalent to retarding the 
density in the original centrally differenced scheme equation (13). This 
leads to 






pWi 


where: 




^i,j.k4“ 




'’i4.j.k ‘’i+r4»j.k 

''i.j.k4 ‘‘i,j,k+t4 


(18a) 


(18b) 

(18c) 

(18d) 


The r, s, and t indices are selected to ensure the ^ , n» and C differencing 
of p in equation (17) are evaluated using upwind formulas. Thus, 

-1 when =• ° 

1 when j I, < 0 

-1 when > 0 

(19) 

1 when < 0 

f-1 > 0 

I 1 when < 0 
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This maintains an upwind influence in the differencing scheme for supersonic 
regions for arbitrary mesh and arbitrary velocity vector. The artificial 
viscosity coefficient v is computed by: 




MAX[(h2^j^^- 1)C, 0] for U.,1 > 0 


( 20 ) 


0] for 


The parameter C is a constant specified by the users. C is usually set 
between 1.0 and 2.0. Use of larger C increases the artificial viscosity 
(upwinding) in the difference scheme. An additional constriant is v £ 1 to 
improve the stability (reference 1). Expressions for v at i,j+l/2,k and 
i,j,k+l/2 can be written similar to equation (20). 


C. AF2 Algorithm 

The AF2 fully implicit approximation factorization scheme (reference 1) for 
three-dimension full -potential equation can be expressed by choosing the N 
operator of equation (2) to be 


‘'^^i.j.k 

- aE;\](« . )CJ^ 


(21a) 


where 


A. 

1 






(21b) 
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(21c) 



The a is a free parameter defined subsequently, and the operator is a 
shifting operator in ^ given by; 

The AF2 scheme can be implemented in a three-step form given by: 


( 22 ) 


Step 1 


(o - 


4 « A.« 

A|^ n J n 






*Vl^i,j,k+l 


(23) 


Step 2 





step 3 




(24) 


(25) 


Here « is a relaxation parameter to be defined by the user with a default 
value of 1.8. Steps 1 and 2 are solved for each k plane in the decreasing k 
direction. This provides f. j ^ for all grid points. Step 3 is then solved 
in the increasing k direction for C. . 

1 9 J fK 

Within each factorized step, the boundary conditions for the intermediate 

solutions (e.g., g and f) should be specified. Initiation of Step 1 at S = 

^max plane reqires information of f at NK+1 which is 

unobtainable. A solution is to set f at NK+1 equal to zero. Similarly, a 

boundary condition is implemented for g at n = (j » 1) and n = 

fmn max 

(j * NO) by imposing (g^)^^! » (9 „)i,nj “ 

For block-structured grid, a boundary condition for f is required at "5 « 
^min ^ “ ^max implemented by 
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imposing j ^ » fmi j ic “ This is consistent with the AF2 scheme 
written in correction form, since convergence implies zero corrections 
everywhere and, therefore, zero f everywhere. 


D. Temporary Damping 

For the N operator in the AF2 scheme, one of the |, n. or S difference 
approximation must be split between two factors. This generates either a 
^ 1 t* ^nt* ^St term, and if it is properly upwind differenced 
(reference 7), provides time-dependent dissipation to the convergence 
process. Following from reference 8 , the 5 difference approximation is split 
between two factors. The control over n (or % ) coordinate direction is 
achieved by adding (or term explicitly in Step 1 (or Step 2). 

Further detail can be found in reference 1. 

The a in equation (21a) can be considered as l/at. This direct analogy to 
time suggests the use of large time step (i.e., small o) to advance time as 
fast as possible for obtaining fast convergence. Reference 9 pointed out that 
this is effective for attacking the low frequency errors but not the high 
frequency errors. An effective approach is to use an a sequence containing 
several values of o. The small values are particularly effective for reducing 
the low frequency errors, and the large values are particularly effective for 
reducing the high frequency errors. A suitable a sequence with analytical 
estimated end points (o|^, o^) is given in reference 9. The sequence is 
generated using the end points to form a geometric sequence. 


E. Boundary Condition 

Along a solid surface, flow tangency boundary condition is imposed. This 
implies that the relevant contravariant velocity component should be zero 
which is simulated by a reflection procedure to account for the flux balance. 
The zeroing of this component also provides a relation to evaluate the 
derivative of d along a grid line normal to the solid surface. This 
eliminates the need to use one-sided difference on d. At the outer boundary 
of the computational mesh, ( J = j = 5 ^ 3 ^. or 5 = 5 the 

d and p are assigned at the freestream values. 
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IV. SOLUTION BY ZONAL DECOMPOSITION 


Block-structured grids are effective for descritizing the computational 
domain. Within each block a grid of good quality can be generated. The 
overall block structuring generally leads to nonanalytical block boundaries 
and nonstandard grid cells. An interesting way to treat this class of grids 
is to solve the flow by zonal decomposition using interfacing grids. The 
problems with the nonanalytical grid block boundaries and nonstandard grid 
cells will be removed with the implementation of suitable inter-block 
continuity conditions between neighboring blocks of grid in the flow solution 
algorithm. In addition, noncontiguous interfacing grids may be used to 
increase the flexibility of the grids. 

A. Model Problem 

The test problem formulated to explore the issues of matching is shown in the 
first part of Figure 6. A rectangular domain AEBCFD was chosen, with the 2D 
Laplace's equation: 

“ ^xx * ^yy ® ® 
defined on a rectangular box 

0£X£2, 0£y£l 

On the box boundary, boundary conditions were chosen of the form 

d = sinh(x)sin(y) + cosh(x)cos(y) (27) 

A boundary EF dividing the box into two domains was defined. Solutions for 
left and right sides were computed independently and matched across EF by the 
procedure discussed below. These solutions were then compared with the exact 
solution calculated for a single domain comprising the large box with the 
boundary EF deleted. 


20 



D 


C 


7 ^^ = 0 
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A E B 


^ given on 
AEBCFD 


Step 1 To start 


FF* 

n — — 

Left box li Right box 

II 
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EE' 


An initial estimation of ^ is assigned along EF and’E'F' where 
^EF “'^E'F'* right box separately. Compute 
(^n^EF ^*^n^E'F' right box. Define 


♦n “ ^&'l*n^EF * ^^n^E'F']* 


Figure 6. Alternating Boundary Conditions Flow Field Matching. 
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Step 2 Neumann Boundary Condition Matching 


D FF' C 

fl 

Left box Right box 

(♦n>EF 

- 


A EE* B 


Neumann Boundary Condition is imposed along EF, (^„)£p * and 
along E*F'» (♦n)£ipi ** ♦p* Again we solve left and right boxes 

separately and define ^E'F'^* 


Step 3 Dirichlet Boundary Condition Matching 


D FF' C 

IT 

Left box II Right box 

•^^EF •Hh^'^E'F' 

II 

W 

A EE' B 


Dirichlet Boundary Condition is imposed along EF and E'F'» 

•^EP ** *^e'F' right boxes are solved separately. 

Define * ^C*^n^EF ^ ^^^n^E'F'^* ^ until 

solution converges. 


Figure 6. (Continued) 
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B. Flow Field Matching Procedure 

The flow field is divided into two boxes as illustrated in Figure 6. The flow 
in each box is solved by an AF scheme. The matching scheme used involved a 
sequential application of alternating boundary conditions (termed an ABC 
method) as outlined in the sequential part of Figure 6. 


C. Calculation of Along EF 

In the first example, grids which 
were contigous at EF were used. For 
either the left or right box, along 
the matching surface, the normal 
derivative need be computed 
using informat ing coming from one 
side of the surface. We note that 
the normal to the matching surface is 
aligned with x. For the left box, 
is computed at node 2 (see 
adjacent sketch). We construct a 
difference approximation to equation 
(26) and require that this equation 
be satisfied at node 2. This leads to 


FF’ 



i (d3-2dp+ dj 
ay ^ 


1 dp~d/i 

O.SaX f^^x^2 ° 


(28) 


solving for { 4^)2 gives 


(dJ 


^2"^4 O.SaX 


x'2 


TT 


ay 


(d3 - 2d2 + dj) 


(29) 
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For the right box at point 6, a similar formula can be derived 


(30) 

When convergence between right and left solutions is achieved, both the 
solutions and its normal derivative should be continuous across the matching 
surface. 


di = dg » ^2 “ ^6 » ^3 “ h 

and (d ) « (d ) 

X 2 X g 

This implies that 

(dg - 2dg * (^3 - 2dg + d^) = 0 


(31) 


This is consistent with the requirement that the finite difference equation 
for Laplace's equation be satisfied at node 2. 


D. AF Scheme 

An iterative scheme for Laplace's equation can be written in correction form 
NC" + uR" = 0 

C" = d"''^ - d" (32) 

r" = Ld" * (« + « )d" 

'XX yy' 

where 6 d and 6 d are the central difference approximations for d^^ 

jj XX 

and dyy respectively, L is the difference approximation of Laplace's 
operator, u is a relaxation factor. The matrix N is factorized such that 
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( 33 ) 


aNC^ . = — (a~6 ) (o~6 )C^ • = — au>Ld^ 

i,j ' yy'^ xx' i,j ^ 

The flow field is solved in two steps. 

Step 1 

(a-6.„)f; ,• = au.L(#; . = oo,r" (34) 

JJ * * J • • J 

The intermediate variable f^ ^ is solved from equation (34). 

• * J 

Step 2 


(o - 6 )C" . 

' xx^ 1,J 



(35) 


The correction vector C. .is solved from equation (35). The m is set equal 

• • J 

to 2. The following formula (reference 10) is used to compute a sequence of 

N a. 


a • 
mm 


= 1 


“max 



. a . n-1 

= (jniii) Nir 

“max 


max 


(36) 


E. Baseline Solutions 

For the boundary condition given by equation 27, the exact solution for 
Laplace's equation is 

^ex * sinh(x) sin(y) + cosh(x) cos(y) (37) 

Numerical solutions for the same problem, on the entire flow field without 
using matching, were also computed on a coarse uniform grid consisting of 10 x 
5 cells and on a fine uniform grid consisting of 20 x 10 cells. 
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For the coarse grid, after 15 iterations, the maximum residual reduced by 
seven orders of magnitude and the maximum error in i (compared with (t^^) was 
0.00174. For the fine grid, it took 16 iterations to achieve the same 
reduction of the maximum residual. The maximum error in d was 0.00045 for 
that case. 

The numerical solution is highly converged and the errors in d reflect the 
truncation error of the finite difference approximation. These numerical 
solutions without using the flow field matching procedure will be referred to 
as baseline solutions. 


F. Flow Field Matching With Contiguous Interfacing Grids 

This simple example is effective for illustrating the approach of flow field 
matching. Refer to Figure 6, a coarse grid with 5 x 5 cells for the left box 
and a coarse grid with 5 x 5 cells for the right box was used. One AF sweep 
was used each time the boundary condition at the matching surface was 
changed. Thus, due to the alternation of boundary conditions, each iteration 
for flow field matching consists of two AF sweeps. It took 24 iterations for 
the maximum residual to reduce by seven orders of magnitude. The maximum 
error in d was 0.00174, which is the same as that in the baseline solution for 
the coarse grid. The same matching problem was also computed on a fine grid 
in which the left and right boxes each had 10 x 10 cells. It took 39 
iterations to reduce the maximum residual by seven orders of magnitude. The 
maximum error in i was 0.00045, the same as that in the baseline solution for 
the fine grid. Comparisons of convergence histories are given in Figure 7 for 
the coarse mesh and in Figure 8 for the fine mesh. Each iteration in the 
matching solutions take twice as many operations as the baseline solutions. 

Thus the basic matching ideas were demonstrated. Fully converged results of 
accuracy equal to that of a single block situation were obtained, but the 
number of iterations required was significantly greater. 


26 





o 20x10 Cells (Single Block) 



Figure 8. Convergence Histories on Fine Mesh. 





6. Flow Field Matching With Noncontiguous Interfacing Grids 

To further explore the versatility of the flow field matching approach, we 
studied a more difficult case where the left box and the right box have 
different grid resolutions. The left box consisted of 5 x 5 cells and the 
right box had 10 x 10 cells (Figure 9). Consequently, the two boxes are 
matched with noncontiguous interfacing grids. This test case is significant 
in that effective application of noncintiguous interfacing grids can reduce 
the number of grid points in the farfield and concurrently reduce the high 
aspect ratio grid in those regions. Grid quality will be improved and total 
number of unknown can be reduced to achieve good numerical accuracy and high 
computational efficiency. For the noncontiguous interfacing grids case, 
special treatment of the flow solution - 
algorithm is needed at the matching 
surface. Refer to Figure 9, along the 
matching surface there are two types of 
nodes. A node which is common to both t\ 
left box and the right box is called a 
common node. A node belonging to the 
right box only is called an uncommon nod< 

For a common node, we first discuss the 
computation of from the solution of 
the Oirichlet boundary condition part of 
the matching problem for point 2 (see 
adjacent sketch). The 4y^ is computed 
using equation 29, where ax and ay are 
replaced by axj and ayj^ respectively. 

For point 6, is computed using 
equation (30) with ax and ay replaced by 
ax 2 and ay 2 « Here 

ax 2 « O.Saxj^ and ay 2 » O.Sayj^ (38) 

The average quantity dj^ used in the iteration procedure outline in Figure 6 
is then obtained through averaging the corresponding values at node 2 and node 
6 . 
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5x5 cells 10 X 10 cells 

for left box for right box 

O Connon node: A node which is cornnon to both the left box 
and the right box. 

□ Uncommon node: A node which belongs to the right box only. 
Figure 9. Flow Field Matching With Noncontiguous Interfacing Grids. 
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The solution of the Neumann boundary 
condition part of the cycle provides (#2 
and tfg, and d Is computed through 
averaging of ^2 and ^ 5 . This defines 
the boundary condition for the Dirichlet 
boundary condition part of the matching 
problem. 


We next consider the uncommon node assuming 
that the solutions of the Dirichlet boundary 
condition problems are given (see adjacent 
sketch). For node 6 , is computed using 
equation (30) while for point 2, d^^ was 
computed by averaging the d^ at node 1 and 
node 3. The can then be defined as the 
average of (d ^^)2 and (dj^)g. 

When d is solved in the Neumann boundary condition part of the cycle, d is 
directly available at node 6 but not at point 2. We may define d 2 by 
averaging d]^ and d 3 > The 7 is then defined by averaging d 2 and dg. 

Numerical experiments show that this scheme encounters a limit cycle. The 
maximum residual can only be reduced by less than two order of magnitude and 
then it starts to oscillate periodically. The oscillation is periodic since 
the sequence of acceleration parameters is changed in a periodic manner. The 
point at which the residual attains its maximum is always at an uncommon 
node. The maximum error in d can only be reduced to the order of 0.01 and 
this error also oscillates periodically. This error also reaches its maximum 
at an uncommon node. 

In order to analyze why the solutions run into a limit cycle, we take a 
further study at the uncommon node. Refer to node 6 , the residual can be 
reduced to zero at node 6 in the right box but not at point 2 in the left 
box. If we force d and to be matched between point 2 and node 6 , then 

there is no guarantee that the residual at node 6 can be reduced to zero in 
the flow field matching procedures. In essence, there are not enough degrees 

I 

of freedom in the system to require both d and dj^ to match exactly at 
uncommon nodes. 
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To correct this problem for the uncommon nodes, we can impose a continuity 
condition for either d or but not both. Numerical experimentation was 
subsequently conducted by matching at the uncommon nodes. Results 
indicate that this correction resolves the limit cycle problem. In addition, 
we used a weighted average procedure for the i (and across the matching 
surface to replace the average procedure. This improved the convergence rate. 

With this correction, the maximum residual reduced by seven orders of 
magnitude in 34 iterations. The residual still reaches its maximum at an 
uncommon node. The maximum error in d is -0.00190 at an uncommon node. This 
is to be expected since at an uncommon node we only enforce the continuity of 
dj^ and we allow d to have a jump across the matching surface. 

Alternatively, one can impose continuity of d at the uncommon nodes and allow 
dj^ to be discontinuous. 


H. Conclusions on Zonal Decomposition 

A summary of the zonal decomposition case study is given in Table 1. 


Table 1 Zonal Decomposition Case Study Summary 


Grid 

Interfacing Grids 
Contiguous 

Convergence in 
Rmax to 10-7 

Error 

One block 10x5 


15 

0.00174 

Two block 5x5+5x5 

Yes 

24 

0.00174 

One block 20x10 


• 16 

0.00045 

Two block 10x10+10x10 

Yes 

39 

0.00045 

Two block 5x5+10x10 

No 

34 

0.00190 
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Our overall conclusion at this point is that the method of matching outlined 
above is feasible. If grids of equal spacing interface at the matching 
boundary, the technique of matching does not cause a loss in accuracy. It 
does require more computer time than a single global solution. 

i 

We have also shown a capability for matching using noncontiguous interfacing 
grids (Figure 9). Our conclusion here is that the local accuracy at the 
matching boundary is of the order associated with the coarser of the 
interfacing grids, which is to be expected. 


I. Remaining Issues 

There remain several issues to explore. One would entail a study of the 
interactions involved when a large number of computational domains interact 
through several matching boundaries. A second issue arises when one uses a 
compressible flow equation. In this case, the suitable matching conditions 
are the continuity of d and the normal mass flux. A third issue arises when 
one considers the use of a governing equation capable of displaying mixed 
character such as the transonic equation. In this case, the appropriate 
boundary conditions governing the solution in a single domain are not known a 
priori, as is the case with the simple Neumann or Dirichlet choices associated 
with Laplace's equation. In the transonic case, the options are: (i) two 

boundary conditions required (when supersonic flow is entering through a 
surface with a direction such that the local Mach cone enters the domain), 

(ii) one boundary condition required (either the local flow is subsonic, or in 
a supersonic case, the local Mach cone cuts the adjacent boundary surface) 
which can be either Neumann or Dirichlet, or (iii) no boundary condition 
allowed (as when the flow is leaving the domain at supersonic speed with no 
portion of the local Mach cone cutting the adjacent boundary surface). Since 
just a prescription of boundary values about a domain does not, in general, 
reveal whether the local direction or Mach number is such that two, one or no 
boundary conditions are allowed, one must begin to think in terms of adaptive 
or selective boundary conditions to be selected during the solution process 
according to the evolution of the local flow. 
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V. GRID EFFECT STUDY ON TWINGE CODE 


Computational experiences of TWINGE reported in reference 1 are based on a 
single-block ring grid. A typical ring grid can be seen in Figure 3. Such a 
grid is very effective for simple geometries, such as a simplified wing. The 
resulting grids for such simple problems are surface adapted, smoothly 
varying, nearly orthogonal, with good cluster control, containing no 
singularities, and the grid cell aspect ratio is moderate everywhere. 

Some of these characteristics may not carry over when we use the 
block-structured grid. Such grid structures can produce fictitious corners, 
non-smooth block boundaries, and grid cells of high aspect ratio. The 
algorithm issues in the flow solution scheme relative to this class of grids 
have not been addressed previously. It was necessary to study these grid 
effects related to TWINGE for adapting the TWINGE code to a block-structured 
grid. 


A. Density Calculation in TWINGE 

In TWINGE, density p is calculated using equation (7d). It follows that; 


r ^ 



V 

^ i 




n 

1 w J 




* h 

* * *2^ * Vc ) ^ 

* (Vt " Vn " > h 

= ''l^l * 

* ) (39) 
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and and to Ag are scaled by the Jacobian determinant D in the 
following manner: 

*2 * "y 
A4 - ( 

Aj-(nyCy)/D .. (40) 

Consistent with the notation of equation (40), the formulas for computing the 

2 

contravariant velocity components U, V, W and q should be written as: 


U 


^ DA^ ^ 



V 

^ s 4 

DA4 A2 DAg 

> 4 

K 



DAg DAg DA3J 


. d^ ^ 


- Aj*^ ^ DCAj*^ ^ * Aj4^4^* Ag<„4^)] 


(41) 


(42) 


1. Approximation Using Averaging Procedures 

In TWINGB, the detailed p computation procedure is the following: 

o- & 0 

i i+1/2 i+1 

d*. , (# » and dr are computed at half i point, i.e., the i+1/2 point shown 

in the adjacent sketch. Aj^ to Ag and D are computed at point i and i+1, 

and then averaged to obtain the corresponding values at i+1/2. Equation (42) 

2 

is then used to compute q at i+1/2. 
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2. Approximation Using Binomial Expansions 
2 

Once q*’ is computed, equation (5d) can be used to compute p. However, in 
TV/INGB equation (7d) is approximated by a five-term binomial expansion. For 
y = 1.4, equation (7d) may be written as 


» - (1 - (43) 

Define 

Q » q^/6 (44) 

p - (1 - 0)2*5 

Binomial expansion of p in terms of Q then leads to 

p = 1 -2.5Q + 1.875 q2 - 0.3125Q^ + 0.0390625Q^ (46) 


Equation 46 offers a much faster p computation than equation (7d). Equation 
(46) is reasonably accurate for engineering application in the transonic 
region. Let p^^ be the exact p computed from equation (43) and p^ be the 
approximation solution of p computed by equation (46), then the relative error 

Err . (., - (47) 

is bounded by 

Err < 10"^ for M < 1.5 (48) 

where M is the local Mach number. 

2 

In order to compute p accurately, we must compute q accurately first. It 
is, therefore, necessary to review the accuracy of the q^ computational 
procedure in TWINGB, focusing on the error introduced by an arbitrary grid. 

Take a uniform flow aligned with the x axis over an orthogonal grid that is 
aligned with the x and y axes. The grid effect can be studied by employing a 
ID model in the x-direction. 
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B. Effect Due to Uneven Grid Spacing 


In TWINGB, for the ID case 
D * 1/x 
U - D 
Aj - 1/x 

. q2 = D Ajd^ 

Let d a X 

then (d^ )i+l/2 “ ^ 

(X^ » (1 + f)/2 

(X{ ),n - f 
(Aj), = 2/(W) 

^^l^i+1/2 “ 

1 + 3f 

“ 2f(l+f)“ ^°M+l/2 
(^^)i+i/2 “ C0.5(l+3f)/(l+f)]^ 

From equation (50) 

2 

d = x, dj^»l, dj^*lis the exact solution. 


From equation (51) f = 1 implies 



Approximate formula recovers the freestream solution » 1 when the mesh 
spacing is uniform. For f # 1, then in general 




* 1 


We plot ®93iost f in Figure 10 to show the effect due to the 

uneven grid on the result. 


We also used the TWINGB code to compute the p corresponding to a uniform flow, 
^ Moo = ‘86, over a flat plat at zero incident angle. Here we know 

that the exact solutions are 


* Q * <1 » P » P everywhere 

A CD 9 

The grid layout for this computation is given in Figure 11. It should be 
noted that along the horizontal lines z = .4, .6 and .8, the grid spacing in x 

jumps across the block boundary. The density computed along these three lines 

is depicted in Figure 12. 

Along these three lines, the p values deviate from the freestream value p « 

CO 

.7083 signficantly near the block boundary. The error in p is resulting from 
the error in the computation on an uneven grid (Figure 12). At all 

other points, the p computed is close to the freestream value. The error in p 

near the block boundary is so severe that some points p is smaller than its 
critical value p* = .634. 

C. Effect Due to Grid Distortion 

Here we want to examine the effect due to grid distortion. The test problem 
we chose here is that of a uniform flow around a circular cylinder. We use a 
smooth grid depicted in Figure 13 to produce the baseline solutions. This 
grid has 60 cells in the circumferential direction with uniform spacing and 10 
cells in the radial direction. In the radial direction, the grid spacing is 
stretched from the cylinder surface to the farfield with a stretching ratio FR 
= 1.2. The farfield boundary is located at 5 radii from the cylinder center. 
The grid in Figure 13 is then distorted by redistributing the grid points 
along a vertical grid line FR = 1.3 for this one line (Figure 14). 
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Figure 10. Grid Effect on Value Conq>uted In TWINGB. 
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Figure 11. A Simple Grid to Test The 6-id Sensitivity of. 7WINGB. 


Figure 12. Error In Density Due to Grid Effect. 
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Figure 14. A Locally Distorted Grid for Circular Cylinder Case* FR*1. 


Comparison of local Mach numbers on the cylinder surface for the baseline grid 
and the distorted grid are given in Figure 15 for freestream Mach number; M^ = 
.45 and in Figure 16 for M,, = .01. The solutions are severely affected near 
the location where the grid is distorted in Figure 15. This effect exists 
even for the case of nearly incompressible flow with M^ « .01 (Figure 16). 


D. Effect Due to Grid Aspect Ratio 

Bloclc-structured grid sometimes results in some grid cells with high aspect 
ratio. It was therefore neccessary to study the sensitivity to grid aspect 
ratio for the basic TWIN6B code. For the cylinder flow, we altered the 
baseline grid of Figure 13 by increasing the grid spacing stretching ratio FR 
to 1.5. 

The resulting grid is given in Figure 17 and numerical experiments produced 
divergent results from TVIIN6B for both M* » .45 and .01. This suggests that 
the stability of the flow solution algorithm in TWIN6B can be severely 
affected by the presence of high aspect ratio grid cells. 

E. Conclusions on Grid Effect Study 

Numerical experiments showed that the inconsistent mapping procedures in 
TWINGB lead to significant error that is proportional to the nonuniformity of 
the grid. The error induced can be severe in the presence of local grid 
distortion. This error can cause instability in the iterative process when 
the grid aspect ratio is high enough. A summary of the TWINGB grid effect 
study can be seen in Table 2. 

This error can be eliminated by using a consistant mapping procedure. The 
difference operators for differencing 4 and for generating the Jacobian matrix 
should be the same to ensure accuracy in the arbitrary grid. E.g., in the ID 
case, all difference formulas should be centered at the Half-I point. 

The extension of this concept to 3D is straightforward if we follow the 
concept of isoparametric mapping taking a suitable subdomain definition (see 
!Appendix A for more detail). Refer to Figure 18, inside the subdomain, 4 and 
X, y, z can be locally approximated by a consistent polynomial approximation. 
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O Smooth Grid 


□ Distorted Grid 


O ® o 


Figure 15. Grid Effect Study. Circular Cylinder. Mo9“.45, Upper Surface M. 




Figure 17. A Smooth Grid for Circular Cylinder Case, FR=1.5. 



Table 2 TWINGB Grid Effect Study Summary 


Test 

Flow Type 

Discussions 

Grid spacing 
discontinuity 

High subsonic 

Discontinuity in the results 

Local Grid 
distortion 

Transonic 

Very poor results locally 

Local Grid 
distortion 

Low subsonic 

Poor results locally 

High Grid 
Aspect Ratio 

Transonic 

Diverged 

High Grid 
Aspect Ratio 

Low subsonic 

Diverged 


Let f be representing either d or x, y, z, then a suitable expansion formula 
for the computational element in Figure 18 may be written as 




a 5 - 0.5 


N 


2=5+ 0.5 


a 0.5 n (n - 1) Mg » 1 - M 3 a 0.5 n (n + 1) 

a 0.5 C ( ;- 1) Lg = 1 - x} L 3 a 0.5 ; { C+ 1) 

and the computational element is defined by 


(52a) 


(52b) 


-0.5 <^<0.5 - 1.0 < n < 1.0 -1.0<^<1.0 (52c) 

The consistent difference formulas for d and x, y, z can be obtained by 
differentiating the polynomial approximations and taking the cell center 
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5 « n « C = 0 as the referencing point for evaluating the difference 
formulas. This improved method can be used in the flow solution algorithm to 
extend TWINGB more effectively to flow problems involving arbitrary grids with 
block structure. 
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VI. . RESULTS FROM ADAPTATION OF TWIN6B TO BLOCK-STRUCTURED GRIDS 

The study of grid effects on TUINGB indicated limitations of TWINGB for 
arbitrary grid. The modified TWINGB (structured for block-structured grid) 
was tested on cases with a typical parabolic arc airfoil and wing using 
three-block grids. These numerical experiments identified the type of 
problems that Can be treated by the current flow solution algorithm in WINGB, 
and suggested that the current algorithm should be improved using a consistant 
mapping procedure (Appendix A). 

The modified TWINGB code was first tested in a 20 mode to solve the airfoil 
problems. The 20 flow was simulated by placing a nonlifting rectangular wing 
between two parallel walls extended to infinity. The wing is non-twisted, and 
consists of same section geometries for all spanwise stations. There are five 
computational stations in the spanwise direction and these stations are 
covered by the same 20 grid. Connecting the corresponding grid points between 
all neighboring stations then produces a 30 grid. The results computed 
between stations were also the same. Only one station is needed to show the 
wing section geometries, grid layout, and surface pressure results for a 
particular test case. 

We then test the modified TWINGB code in a 30 mode on wing problems. For this 
study, the rectangular wing planform in the 20 case was changed by adding 
sweep to the planform. This generates a swept non-tapered wing between two 
parallel walls extended to infinity. To study the 30 effects, including the 
variations of the solution in the spanwise direction, 20 computational 
stations were used in this direction. For simplicity, each station is still 
covered by the same 20 grid. 

Unless otherwise stated, the parameters for AF2 iterations are set at default 
values. In the default mode, the time-like dissipation coefficient p,. .1, 

e « 0, the artificial viscosity coefficient C =1.2, the relaxation factor 
u = 1.8, the a sequence end points are (oj_ = 0.4, oj^ = 4), and there are 
eight elements in the o sequence. 
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All TWINGB calculations for the 20 cases use a mesh with 6,000 cells (80 x 5 x 
15). This places 40 cells from leading to trailing edge. The TWINGB 
calculations for the 3D cases use a mesh with 22,800 cells (80 x 19 x 15). 

The grid spacing in the spanWise direction (y direction) is uniform for both 
20 and 30 cases. 


A. Test Cases for Airfoil Problems 

1. Uniform Flow Over a Flat Plate 

For this test case, the grid in each spanwise station has 80 cells in the x 
direction (streamwise direction) and 15 cells in the z direction (normal to 
the wing direction). For a grid with equal spacing in x and equal spacing in 
z (Figure 19), numerical experiments at = .9 showed that the AF2 iteration 
is stable and that it completely recovers the freestream flow. Only the lower 
half space is needed since the flow is symmetric with respect to z = 0. The 
lower half space (instead of the upper half space) should be used to make 
the 4 » n, and C coordinate system a right hand system - the same system as 
the X, y, and z coordinate system. The use of inconsistent coordinate systems 
would lead to numerical instability in the AF2 iteration beyond the simple 
case of two dimensional totally subsonic flow. 

We next stretched the grid spacing in the z direction with a stretching ratio 
FZ s 1.2. The resulting grid is depicted in Figure 20. This grid leads to 
numerical instability with TWINGB for both = .9 and .0001. The maximum 
residual and maximum correction versus iteration is plotted in 

Figure 21 for the case of = .0001. Clearly, the iteration is diverging. 

The peaks in the AF2 convergence (or divergence) history correspond to the 
small value of o (i.e., larger pseudo-time step) in the eight element sequence. 

We then reduced FZ to 1.1 (Figure 22), and the resulting AF2 iteration for the 
flow solution becomes stable (for the case of = .9) and the solution 
recovers the freestream flow. 
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Figure 19. 6rid-A. Flat Plate, Uniform Grid. 
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Figure 21. Growth of and for Flow Over a Flat Plate, Grid-B, 
M»*.0001, 
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We next generated a three-block grid keeping FZ=1.1 (Figure 23). The first 
block (left hand block) contains i * 1 to 21, the second block (central block) 
i s 21 to 61, and the third block (right hand block) i = 61 to 81. For the 

central block, uniform grid spacing in x is used. For the left hand and right 

hand block, the x grid spacing was stretched with a stretching ratio FX = 
1.125. The block boundaries are located at x « 0 (i = 1), x x= .4056 (i « 21), 
x = .5944 (i » 61), and x = 1 (i =81). The flat plate is placed at z = 0 

from X = .4056 to x = .5944 with a chord length of .1888. The far field in 

the 2 direction is at z = -1. With this grid, the AF2 iteration is stable and 
the convergence histories of and are given in Figure 24 for 
= .9 and in Figure 25 for M = .0001. 

The exact solution for this case is Cp =0 everywhere. To check the 
solution accuracy based on the grid in Figure 23, the pressure coefficient 
Cp on the flat plate surface is plotted against the normalized x coordinates 
(Figure 26) for the case of M^= .90. The resulting Cp is very small but 
not zero. This illustrated the grid effects on the AF2 flow solution 
algorithm. This test case is simple and the relevant results appear to be 
acceptable within engineering accuracy. It remains to be examined what grid 
effects will do for more difficult cases. 


2. 0.03 Thick Parabolic Arc Airfoil 

The three-block grid of Figure 23 can be perturbed slightly to produce a grid 
for the case of flow over a thin parabolic arc airfoil. A grid so produced 
(Grid-E) is depicted in Figure 27 for a 0.03 thick parabolic arc airfoil. The 
freestream Mach number M^^ = .75 and the convergence history for and 
are shown in Figure 28. Rapid reduction of and are 
realized showing the effectiveness of the AF2 iteration for this test case. 

The which is a measure of how the governing equation is satisfied, 

drops three order of magnitudes in 37 iterations. 

The surface Cp is shown in Figure 29 (square symbol). In the same figure, 
the Cp computed from TCAS is given by the solid line for comparison. The 
TCAS code was modified from TAIR code by Kwak (private communications at NASA 
Ames). Both codes produce symmetrical results - a correct flow behavior for 
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Iterations 


Figure 24. Convergence Histories for Flow Over a Flat Plate, Grid-D, 
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Figure 25. Convergence Histories for Flow Over a Flat Plate, Grid-D, 
M*«.0001. 
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2D subcritical flow over a symmetrical airfoil. The general agreement between 
the two results is good. The differences in the two solutions are mainly 
caused by two factors. Firstly, different grids are used; specifically, TCAS 
uses many more grid cells in the z direction. Secondly, some detail in the 
solution procedures are different between the two codes (e.g., p is computed 
at cell centers for TCAS but computed at half-I points for TWINGB). 


3. 0.12 Thick Parabolic Arc Airfoil 

The grid of Figure 27 can be further perturbed to produce a grid for flows 
over a 0.12 thickness parabolic arc airfoil (Grid-F, Figure 30). Grid-F was 
used for calculating the transonic case of M_ « .8. 

Increasing the airfoil thickness produces stronger perturbation to the uniform 
flow. The flow becomes transonic for this combination of airfoil thickness 
and freestream Mach number. The solution process for this case was divergent 
The growth of and are given in Figure 31 for * .8. The AF2 
iteration appears to be slowly divergent in an oscillating manner. Adjustment 
of parameters in the AF2 iteration seems to be ineffective in getting the 
solution to converge. 

The growth of and are given in Figures 32 to 35 to study the 
parameter variations. Figure 32 is the case for » = 1, apparently the 
solution is not convergent. We then increased the magnitude of the a sequence 
end points to = 1 and = 20. This has the effect of reducing the 
time-like step to improve the numerical stability. Even with such an a 
sequence, the solution is still not convergent. We next increased the 
time-like dissipation coefficient 6^ from .1 to 1 and add it everywhere in 
the flow field. The solution is still not convergent (Figure 34). The 
numerical instability appears to be caused by two major effects, the transonic 
effect and the grid effect. In order to solve this case, we need to study 
each effect separately. 

These effects were studied individually as follows. Using the Grid-F for the 
case of = .0001, the convergence histories are given in Figure 35 where 
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Figure 27. Grid-E, .03 Parabolic Arc Airfoil, FX=1.125, FZ=1.1. 




Figure 29. Surface Cp for Flow Over a .03 Parabolic Arc Airfoil, 
6rid-E. H -.75. 
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Figure 31. Growth of and for Flow Over a ,12 Parabolic Arc 
Airfoil, 6r1d-F, 
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the solution is converging. Apparently the transonic calcuations are more 
sensitive to the error induced by the grid effect. To make the transonic case 
work, it is clear that the 6rid-F needs to be improved. Examining the central 
block (from i » 21 to i » 61)' to the grid in 6rid-F, one sees that the grid 
aspect ratio near the far field (z * -1) is quite large. 

A method to reduce the grid aspect ratio near z » -1 is to reduce FZ, e.g., FZ 
» 1 so as to make the grid spacing in z uniform. In so doing, the grid aspect 
ratio near the airfoil surface is getting worse (6rid-G, Figure 36). This 
grid is not good enough to make the solution converge for » .8 (Figure 
37). A better method to reduce the grid aspect ratio at z » -1 is to 
redistribute the grid points along z = -1 (e.g., uniform grid spacing in x 
along z = -1, grid-H, Figure 38). 

For this grid, we also readjusted the block boundaries slightly so that the 
airfoil leading edge is at x « .4, and trailing edge at x = .6. The 
convergence histories in Figure 39 shows that the is reduced by three 
order of magnitude in 119 iterations. In Figure 39, R_j,„ and are 
plotted for every 8 iterations, corresponding to a « 0 |_ in the eight-element 
sequence. Also shown in Figure 39 is the convergence history for the number 
of supersonic points NSUP. The NSUP tends to be overshoot within the first 
eight a elements of the iterations, but it converges to its final value in 
only 27 iterations. This shows that the size of the supersonic pocket is 
stabilized very fast. This is typical for AF2 scheme in transonic flow 
calculations 

The surface Cp distribution is given in Figure 40 and compared with the 
solutions using TCAS and LTRAN2. LTRAN2 used a small disturbance 
formulation. This is the main cause of the deviation of the LTRAN2 solution 
from the other two solutions. General agreement is quite good and the 
magnitude of the small difference between the solutions of TCAS and TWINGB is 
the same as for the subcritical case. 


B. Test Cases For Wing Problems 
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Figure 36. Grid-G. .12 Parabolic Arc Airfoil. FX=1.125, FZ=1. 
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Figure 37. Growth of and for Flow Over a .12 Parabolic Arc 
Airfoil, 6r1d-6. Mj^.8. 
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Figure 38. Grid-H, .12 Parabolic Arc Airfoil. FXl-1.129. FX2-1, FZ-1.1. 
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Figure 39b. 


Convergence History of NSUP for Flow Over a .12 Parabolic Arc 
Airfoil, 6rid-H. M^.8. 
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Figure 40. Surface Cp Distributions for Flow Over a .12 Parabolic Arc 
Airfoil. Grid-H. M.-.8. 
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1, Uniform Flow Over a Swept Flat Plate. 

This case was not tested in view of the success of the cases of flow over 
parabolic arc airfoils. 

2. A Swept Wing With a .03 Parabolic Arc Wing Section. 

A high aspect ration (AR = 9.5) wing has. been studied first. This wing is 
non-tapered, has a 30 degree sweep and is mounted between two parallel walls 
extended to infinity. The planform of this wing is identified as planform-A. 

In each of the 20 spanwise stations, a grid similar to that in Figure 27 
(6rid-E) is used. The block boundaries are adjusted slightly. For the wing 
root section, the wing leading edge is at x = .4 and trailing edge at x = .6. 

This wing is analyzed at = .75. The convergence histories given in Figure 
41 shows that it takes 103 iterations to reduce the by three orders of 
magnitude. The 30 effects slow down the convergence and the resulting 
convergence speed is two and a half times slower than the 20 counterpart. 

This is to be expected since the variation of the flow in the spanwise 
direction may induce further error due to the approximations employed in 
TWIN6B. Since these errors are grid dependent, the convergence should be 
improved when using a better grid. An improved grid (compared with 6rid-E) is 
displayed in Figure 42 as Grid-I. Using Grid-I for flow calculations, the 
convergence histories given in Figure 43 shows that reduced by three 
orders of magnitude in only 54 iterations. The convergence speed is improved 
by a factor greater than two due to grid improvement. The surface Cp 
results are given in Figure 44 for the wing root, wing mid-span, and wing tip 
stations. 

The low aspect ratio (AR = 1.9) case has also been analyzed at = .75 and 
this case is expected to produce a more pronounced 3D effect. The planform of 
this wing is called planform-B which has a 30 degree sweep. The wing is 
non-tapered and it is mounted between two parallel walls extended to 
infinity. In each of the 20 spanwise stations, Grid-I (Figure 42) is used. 
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Figure 41. Convergence History for Flow Over a Swept Wing With a .03 
Parabolic Arc Wing Section, Planform-A, Grid-E, M„=.75. 
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The numerical convergence for this case is very slow (274 iterations, see 
Figure 45). The pronounced 3D effects affect the convergence and the 
resulting convergence speed is two arid a half times slower than the high 
aspect ratio wing case. The surface Cp distributions are given in Figure 46 
for the wing root, wing mid-span, and wing tip stations. 

3. A Swept Wing with a 0.12 Thick Parabolic Arc Wing Section 

For the high aspect ratio wing case, planform-A (AR = 9.5) is used. In each 
of the 20 spanwise stations, 6rid-H (Figure 38) is used. 

The wing is analyzed at » .8. The convergence histories given in Figure 
47 show that it takes 123 iterations for the to reduce by three order 
to magnitudes. However, the number of supersonic point NSUP converged in 53 
iterations. This shows that the size of the supersonic pocket is established 
at very early stage of the iterations, similar to the 2D cases. 

The shock wave position is plotted on the wing planform in Figure 48. As 
expected, the shock wave approaches both sidewalls nearly perpendicularly. 

Near center span, the shock is swept and is approximately parallel to the wing 
leading and trailing edges. The surface Cp results in the root, center span 
and tip planes are given in Figure 49. 

For the low aspect ratio wing case, planform-B (AR = 1.9) is used. Grid-H is 
used for each of the 20 spanwise stations. The wing is analyzed at « .8. 

Numerical convergence is extremely slow (Figure 50). The convergence of 
Rj^ax three orders of magnitude takes 373 iterations and the convergence 
of NSUP consumes 135 iterations. That is more than twice as slow as the high 
aspect ratio wing case. This is resulting from the error, induced from the 
approximation procedures in the flow solution algorithm, becoming more severe 
due to strong 3D effects. 

The shock wave position is plotted on the wing planform in Figure 51. The 
shock wave approaches both sidewalls nearly perpendicularly. Near center 
span, the shock is swept and is approximately parallel to the wing leading 
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Figure 47a. Convergence History of for Flow Over a Swept Wing With a 

.12 Parabolic Arc Wing Section, Planform-A, Grid-H, AR*9.5. 
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Figure 50b. Convergence History of NSUP for Flow Over a Swept Ming With a 

.12 Parabolic Arc Wing Section, Planform-B, 6rid-H, M^*.8, AR»1.9. 
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and trailing edge. This is qualitatively the same as the high aspect ratio 
wing case. The surface Cp distribution in the root, center span and tip 
plans are given in Figure 52. 


C. Test Case Study Summary 

Results of this chapter fully back up the grid effect study on TWINGB code 
given in Chapter IV where a method to improve the TWINGB flow solution 
algorithm is also presented. It is neccessary that this improved method be 
adapted for extending TWINGB more effectively to flow problems involving 
arbitrary grids with block structures. 

A summary of TVlINGB test case study is given in Table 3. 
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Table 3 TWIN6B Case Study Summary 


Airfoil Problem 

Thickness 

M 

Grid 

Convergence in 

Convergence 

Ratio 



*^max ^0 

in NSUP 

0.03 

0.75 

Grid-E 

37 


0.12 

0.80 

Grid-F 

Di verge 


0.12 

0.0001 

Grid-F 

Converging, terminated 





at 50 iterations 


0.12 

0.80 

Grid-G 

Diverge 


0.12 

0.80 

Grid-H 

119 

27 


30* Swept Wing ARs9.5 


Thickness 

M 

Grid 

Convergence 

Convergence 

Ratio 



^max 

in NSUP 

0.03 

0.75 

Grid-E 

103 


0.03 

0.75 

Grid-I 

54 


0.12 

0.80 

Grid-H 

123 

53 


30* Swept Wing AR«1.9 


Thickness 

M 

Grid 

Convergence 

Convergence 

Ratio 



^ax 

in NSUP 

0.03 

0.75 

Grid-I 

274 


0.12 

0.80 

Grid-H 

373 

135 
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VII. CONCLUSIONS 


This study explores means for extending AF schemes to the treatment of 
transonic flows requiring the use of complex, multi-block grids. A method 
involving "zonal decomposition" with "interfacing grids" was developed. Test 
case studies for imcompressible flow problems are completed and successful. 
This method is judged to be feasible for extension to transonic flows wherein 
the block boundaries are everywhere subcritical. The case of mixed flow along 
a block boundary will be much more difficult.. 

This study also adapted a fully implicit AF based 3D transonic code - called 
TVlINGB to a block structured grid. Subsequent numerical experiments on flow 
over parabolic arc airfoils and wings using block-structured grids showed that 
the flow solution algorithm in TWINGB is sensitive to grid effects. The 
present conclusion is that the usefulness of TWINGB in the present form is 
limited to problems having only high quality grids. Uneven grid spacing or 
high aspect ratio grid cells lead to difficulties. An improved method to 
reduce the grid dependency of the flow solution algorithm in TWINGB is 
presented. The analysis and test case results in this study suggest that it 
is neccessary to adapt the improved method for extending TWINGB more 
effectively to flow problem involving arbitrary grids with block structure. 
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VIII. RECOMMENDATIONS , 


The experiences of this study results in the following four recommendations: 

A. Compatibility between flow solution algorithm and computational grids 
must be carefully studied. This has far-reaching importance when developing 
computational transonic aerodynamic methods for flow over realistic airplane 
configurations using arbitrary grids. 

B. Many approximation procedures in TWINGB are acceptable only for grids 
with very good quality (e.g., globally smoothly varying, moderate grid aspect 
ratio). To extend TWIN6B to complex configurations with more general grids, 
it is necessary that the approximation prpcedures in TWINGB be improved. 
Specifically, the mapping procedure between the physical space and the 
computational space must be fully consistent. 

C. Block-structured or segmented grids offer a sensible way to 
descritize the computational domain. Within each block, a grid of very good 
quality can be generated. The assemblege of these blocks provides the 
complete grid for computation. Generally, this leads to nonanalytical grid 
block boundaries and nonstandard grid cells. A logical treatment for this 
class of grids is to solve the flow by zonal decomposition using interfacing 
grids. In so doing, the problems with the nonanalytical grid block boundaries 
and nonstandard grid cells are circumvented by imposing suitable inter-block 
continuity conditions between neighboring blocks of grid in the flow solution 
algorithm. The success of applying the zonal decomposition method on the 
model problem in Chapter IV warrants a further study of this method for 
transonic wing problems. 

D. Refinement of the multi-block grid generation is recommended. 
Specifically, linear equations for block grid generation should be replaced by 
nonlinear ones. This would ease the control of interior grid points in each 
block and therefore produce a grid with a desirable inter-block continuity of 
the grid points. The effect due to this improvement is expected to be 
significant when block boundaries are highly curved. 


98 



IX. REFERENCES 


1. Holst, T. L., 1979, !'A Fast, Conservative Algorithm for Solving the 
Transonic Full-Potential Equation", AIAA , Vol. 18, 1980, pp. 1431-1439. 

2. Lee, K. D., "3-D Transonic Flow Computations Using Grid Systems with Block 
Structure", AIAA, Fifth Computational Fluid Dynamic Conference, Palo Alto, 
California, June 22-23, 1981. 

3. Lee, K. D., and Rubbert, P. E., "Transonic Flow Computations Using Grid 
Systems with Block Structure", Proceedings , the 7th International 
Conference on Numerical Methods In Fluid Dynamics, Stanford University, 
California, June, 1980. 

4. Thompson, J. F., Thames, F. C., and Mastin, C. M., "Automatic Numerical 
Generation of Body-Fitted Curvilinear Coordinate System for Field 
Containing Any Number of Arbitrary Two Dimensional Bodies", Journal of 
Computational Physics ,, Vol. 15, 1974, pp 299-319. 

5. Yu, N. J. , "Grid Generation and Transonic Flow Calculations for 
Three-Dimensional Configurations", AIAA Paper 80-1391, AIAA 13th Fluid and 
Plasma Dynamics Conference, Snowmass, Colorado, July, 1980. 

6. Middlecoff, J. F. and Thomas, P. D., "Direct Control of the Grid Point 
Distribution in Meshes Generated by Elliptic Equations", AIAA Paper 
79-1462, 1979. 

7. Jameson, A., "Iterative Solution of Transonic Flows Over Airfoils and 
Wings, Including Flows at Mach 1", Communications on Pure and Applied 
Mathematics , Vol. 27, 1974, pp 283-3D9T 

8. Holst, T. L., "An Implicit Algorithm for the Conservative Transonic Full 
Potential Equation Using an Arbitrary Mesh", AIAA Paper 78-1113, July, 
1978. 

9. Holst, T. L. and Ballhaus, W. F., "Fast Conservative Schemes for the Full 

Potential Equation Applies to Transonic Flows", AIAA, Vol. 17, 1979, pp. 
145-152. 

10. Ballhaus, W.F., Jameson, A., and Albert, T. , "Implicit Approximate- 
Factorization Schemes for the Efficient Solution of Steady Transonic Flow 
Problems", AIAA J., Vol. 16, 1978, pp. 573-579. 


99 



APPENDIX A 


CONSISTANT MAPPING AND VELOCITY CALCULATION 

Theorem. Given a velocity potential for a uniform flow, in order to produce a 
uniform velocity, it is necessary and sufficient that the mapping procedure be 
fully consistent. 

Necessary Proof 

Let X, y, z be the physical coordinates and 5 , n» $ be the computational 
coordinates and the mapping between the physical and computational spaces is 
one-to-one, we have 



In equation (Al), the 5* n» and 5 differentiations of d, x, y, and z are 
usually approximated by difference formulas in the n, 5 coordinates. We 
thus approxiamte equation (Al) by equation (A2) 



where subscript d denotes that the d derivatives in |, n and S are 
approximated by differences formulas for d, and subscript xyz denotes that the 
X, y, and z derivatives in |, n and $ be approximated by difference 
formulas for x, y, and z. 

If d = X , then 



(A3) 


9 


d s d s 0 
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and 



Also, from equations (A2) and (A3), if we indeed produce a uniform flow, then 



From equations (A4) and (A5) 



Thus the difference formula for i and x must be consistent. Similarly, we 
can set d * y and d » z to show that the difference formulas for d and x, y, z 
must be consistent. 

Sufficient Proof 
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Assuming the difference formulas for d and x, y, z are inconsistent, we may 
replace the subscripts d and xyz by a subscript D, We then rewrite equation 
(A2) as 


(A7) 


if d » X , then 



If the mapping matrix is nonsingular, then equation (A8) has a unique solution 




This is the uniform flow sought. 


(A9) 
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